Zbadać częstotliwość oraz skalę opóźnień pociągów w danym okresie oraz zależność ich od sytuacji meteorologicznej. Przedstawić tę zależność za pomocą wizualizacji.
Przed badaniem robimy przypuszczenie, iż im gorsza sytuacja pogodowa, tzn. im więcej opadów pewnego dnia oraz im mniejsza temperatura, tym opóźnienia będą większe i częstsze.
Opis: The dataset contains information about delays on polish railroads published by the Polish State Railways in real time. Data was scraped from https://infopasazer.intercity.pl. It covers the period from the midnight of 2022-05-16 until the midnight of 2022-05-30. Observations were collected every 5 minutes.
Zmienne:
datetime - timestamp at which the sample was collected (Warsaw time)
id - train ID
carrier - the name of the carrier (mainly its PKP Intercity)
date - date of the train departure
connection - beginning and the destination for the train
arrival - planned arrival time at the destination
delay - current estimated delay
name - train station name
Dataset został opublikowany w roku 2022 i do teraz ma nie więcej niż 800 pobrań. Nie znalazłem również żadnych raportów z wykorzystaniem tego zbioru, więc można stwierdzić, że zbiór jest dość unikatowy.
Z drugiej strony, zbiór liczy około 4 mln obserwacji, więc jest dla nas dość reprezentatywny.
Opis: Są to archiwowe dane pogodowe dla okresu, odpowiadającego okresu zbierania danych w pierwszym zbiorze. Za lokalizaję wzięto przybliżone koordynaty geograficznego centrum Polski.
W tym projekcie po czyszczeniu zbioru danych oraz konwertowaniu poszczególnych zmiennych do właściwego dla nas typu zajmiemy się badaniem zależności pomiędzy pewnymi typami zmiennych za pomocą np. wyznaczenia korelacji pomiędzy zmiennymi numerycznymi. Znajdziemy też średnie opóźnień dla poszczególnych przewoźników lub relacji, uwzględnimy dane pogodowe, żeby wyłapać tę zależność oraz wesprzymy nasze wnioski odpowiednimi wizualicjami. Pomysłem też jest zrobienie liniowego modelu predykcyjnego, który na podstawie danych pogodowych będzie przewidywał średnie opóźnienie na poszczególnych stacji.
Zacznijmy od importowania zbiorów oraz wszystkich potrzebnych bibliotek
#rm(list = ls())
dftemp = read.csv("delays.csv")
lockBinding("dftemp", globalenv())
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(chron)
##
## Attaching package: 'chron'
##
## The following objects are masked from 'package:lubridate':
##
## days, hours, minutes, seconds, years
library(ggplot2)
library(fuzzyjoin)
library(corrplot)
## corrplot 0.95 loaded
weather = read.csv("open-meteo-52.55N13.41E38m.csv")
df = dftemp
Zajmijmy się najpierw zbiorem opóźnień
summary(df)
## datetime id carrier date
## Length:3718170 Length:3718170 Length:3718170 Length:3718170
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## connection arrival delay name
## Length:3718170 Length:3718170 Length:3718170 Length:3718170
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
sum(is.na(df))
## [1] 0
Konwertujemy typy danych
df$datetime = strptime(df$datetime, format = "%Y-%m-%d %H:%M")
df$date = strptime(df$date, format = '%Y-%m-%d')
df$arrival = strptime(df$arrival, format = '%H:%M')
split_data = do.call(rbind, strsplit(df$connection, " - "))
df$start = split_data[,1]
df$end = split_data[,2]
rm(split_data)
df = df[, -5]
df$id = as.factor(df$id)
df$carrier = as.factor(df$carrier)
df$delay = as.numeric(gsub(" min", "", df$delay))
df$name = as.factor(df$name)
df$start = as.factor(df$start)
df$end = as.factor(df$end)
Teraz czyscimy zbior danych pogodowych
head(weather)
## latitude longitude elevation
## 1 52.54833 13.407822 38.0
## 2 time temperature_2m (°C) relative_humidity_2m (%)
## 3 2022-05-16T00:00 12.8 69
## 4 2022-05-16T01:00 12.0 70
## 5 2022-05-16T02:00 11.6 72
## 6 2022-05-16T03:00 10.5 76
## utc_offset_seconds timezone timezone_abbreviation
## 1 7200 Europe/Berlin CEST
## 2 rain (mm) surface_pressure (hPa) wind_speed_10m (m/s)
## 3 0.00 1016.7 4.07
## 4 0.00 1016.7 3.79
## 5 0.00 1016.5 3.04
## 6 0.00 1016.5 3.04
Widzimy, że jest niedobrze, pierwsze 2 usunąć
colnames(weather) = c(weather[2, 1], weather[2, 2], weather[2, 3], weather[2, 4], weather[2, 5],weather[2, 6])
weather = weather[c(-1, -2), ]
sum(is.na(weather))
## [1] 0
weather$time = strptime(weather$time, format = '%Y-%m-%d T %H:%M')
weather$`temperature_2m (°C)` = as.numeric(weather$`temperature_2m (°C)`)
weather$`relative_humidity_2m (%)` = as.numeric(weather$`relative_humidity_2m (%)`)
weather$`rain (mm)` = as.numeric(weather$`rain (mm)`)
weather$`surface_pressure (hPa)` = as.numeric(weather$`surface_pressure (hPa)`)
weather$`wind_speed_10m (m/s)` = as.numeric(weather$`wind_speed_10m (m/s)`)
NaNów nie ma, dobrze
Lączymy zbiory
df = fuzzyjoin::difference_left_join(
df, weather,
by = c("datetime" = "time"), # Adjusting column names here
max_dist = as.difftime(30, units = "mins") # Adjust max_dist as needed
) %>%
group_by(datetime) %>%
filter(abs(difftime(datetime, time, units = "secs")) == min(abs(difftime(datetime, time, units = "secs")))) %>%
ungroup() %>%
select(-time) # Remove the redundant time column from weather
rm(weather)
Robimy delay factor
df$delay_factor = cut(df$delay,
breaks = c(-Inf, 0, 5, 30, 60, Inf),
labels = c("0 min", "1-5 min", "6-30 min", "31-60 min", "> 60 min"),
right = TRUE)
delay_over_60_by_carrier = df[df$delay_factor == "> 60 min",] %>%
group_by(carrier)%>%
summarise("Delays_Over_60_mins" = length(delay_factor))
library(ggridges)
library(plotly)
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
delay_summary = df %>%
group_by(carrier, delay_factor) %>%
summarise(count = log(n()), .groups = 'drop')
total_counts = delay_summary %>%
group_by(carrier) %>%
summarise(total_count = sum(count), .groups = 'drop')
delay_summary = delay_summary %>%
left_join(total_counts, by = "carrier") %>%
mutate(count = 100*count / total_count)
bar_plot = ggplot(delay_summary, aes(x = carrier, y = count, fill = delay_factor)) +
geom_bar(stat = "identity", position = "dodge") +
labs(
title = "Counts of Train Delays by Carrier",
#x = "Carrier",
y = "Log(Percentage of delay type)",
fill = "Delay Factor"
) +
theme_minimal(base_size = 15) +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
axis.text.x = element_text(face = "bold"),
legend.position = "right"
)
interactive_plot = ggplotly(bar_plot)
interactive_plot = interactive_plot %>%
layout(
xaxis = list(tickangle = 45)
)
interactive_plot
df_by_hour = df%>%
group_by(format(datetime, "%Y-%m-%d %H"))%>%
summarize(mean(delay), mean(`rain (mm)`), mean(`temperature_2m (°C)`), mean(`relative_humidity_2m (%)`), mean(`wind_speed_10m (m/s)`))
df_by_hour = df_by_hour[,-1]
correlation = cor(df_by_hour)
corrplot(correlation)